imjviuufimi 


% 


reCORt^/KTlON  PAGE 

AD-A214  490  k]  irVVUj.  0  "BT" 


Form  Approved 
OMB  Ho  0704-01 8B 


1b  RESTRICTIVE  MARKINGS 


2b.  DECLASSIFICATION /DOWNGRADING 


4  PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 


^LECTEWt 

NOV  20  B  DISTRIBUTION /AVAILABILITY  OF  REPORT 

®  "Approved  for  public  release; 
distribution  unlimited. 


6b.  OFFICE  SYMBOL 
(If  applicable) 


5.  MONITORING  ORGANIZATION  REPORT  NUM8ER(S) 

.  'Jis  r  \ ■ 


7a.  NAME  OF  MONITORING  ORGANIZATION 


6c.  ADDRESS  (City,  State,  and  ZIP  Code) 

Department  of  Electrical  Engineering 
Information  Systems  Laboratory 
rd,  CA  94305-*055 


8a.  NAME  OF  FUNDING /SPONSORING 
ORGANIZATION 


8b.  OFFICE  SYMBOL 
(If  applicable) 


AFOSR  _ 


8‘c  ADDRESS  (City,  State,  and  ZIP  Code) 

Building  410 


Air  Force  _Qf 


7b  ADDRESS  (City,  State,  and  ZIP  Code) 
Building  410 

Bolling  AFB ,  DC  20332-6448 


9  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUM8ER 
AFOSR-88-0327 


10  SOURCE  OF  FUNDING  NUM8ERS 


DC  20332-6448 

PROGRAM 
ELEMENT  NO. 

PROJECT 

NO. 

TASK 

NO 

• 

61102F 

2304 

A6 

WORK  UNIT 
ACCESSION  NO 


1 1 .  TITLE  (Include  Security  Classification) 

STUDIES  IN  STATISTICAL  SIGNAL  PROCESSING 


12.  PERSONAL  AUTHOR(S) 
Professor  Thomas  Kailath 


14.  DATE  OF  REPORT  (Year,  Month,  Day)  115.  PAGE  COUNT 


COSATI  CODES 


GROUP 


SUB-GROUP 


18.  SUBJECT  TERMS  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 


1 9.  ABSTRACT  ( Continue  on  reverse  if  necessary  and  identify  by  block  number) 

This  is  an  annual  report  on  Contract  AFOSR-88-0327,  for  the  period 
July  1,  1988  to  June  30,  1989.  Section  I  provides  a  brief  overview 
of  our  work,  while  the  remaining  sections  describe  in  some  detail 
our  recent  results  on  efficient  factorization  of  structured  matrices, 
polynomial  zero- locat ion ,  and  recursive  layer  peeling.  the  appendices 
list  publications  and  other  activities  during  the  last  year.  The 
abstract  of  the  recently  completed  Ph.D  dissertation  of  J.  Chun  is 
also  appended  to  this  report. 


20.  DISTRIBUTION /AVAILABILITY  OF  ABSTRACT  21  ABSTRACT  SECURITY  CLASSIFICATION 

□  UNCLASSIFIED/UNLIMITED  (g  SAME  AS  RPT  □  0TIC  USERS  UNCLASSIFIED 


22a.  NAME  OF  RESPONSIBLE  INOIVIOUAL  22b  TELEPHONE  (Include  Area  Code)  22 c  OFFICE  SYMBOL 

DR.  EYTAN  BAROUCH _ _ _ (202)  767-4940  NM 


OD  Form  1473,  JUN  86  Previous  editions  are  obsolete.  _ secl 


*  1 1  r.  * 


O 


II 


S  PACE 


i  OSR-TK.  8  9-1353 


AFOSR88-0327  (Task  NO.  2304/A6) 


STUDIES  IN  STATISTICAL  SIGNAL  PROCESSING 


Annual  Report 

Professor  Thomas  Kailath 
Information  Systems  Laboratory 
Department  of  Flectrical  Engineering 
Stanford  University 
Stanford,  CA  94305-4055 


August  1,  1989 


Period  Covered:  July  1,  1988  -  June  30,  1989 


Approved  for  Public  Release;  Distribution  Unlimited 


Ptepared  for: 

Air  Force  Office  of  Scientific  Research 
Bldg.  410 

Bolling  Air  Force  Base 
Washington,  D.C.  20332-6448 


STUDIES  IN  STATISTICAL  SIGNAL  PROCESSING 


This  is  an  annual  report  on  Contract  AF-88-0327,  for  the  period  July  1,  1988  to 
June  30,  1989.  Section  I  provides  a  brief  overview  of  our  work,  while  the  remaining 
sections  describe  in  some  detail  our  recent  results  on  efficient  factorization  of 
structured  matrices,  polynomial  zero-location,  and  recursive  layer  peeling.  The 
appendices  list  publications  and  other  activities  during  the  last  year.  The  abstract  of  the 
recently  completed  Ph.D.  dissertation  of  J.  Chun  is  also  appended  to  this  report. 

1.  INTRODUCTION 

^  The  primary  objective  of  our  research  is  to  develop  efficient  and  numerically 
stable  algorithms  for  nonstationary  signal  processing  problems  by  understanding  and 
exploiting  special  structures,  both  deterministic  and  stochastic,  in  the  problems.  We 
also  strive  to  establish  and  broaden  links  with  related  disciplines,  such  as  cascade  filter 
synthesis,  scattering  theory,  numerical  linear  algebra,  and  mathematical  operator  theory 
for  the  purpose  of  cross  fertilization  of  ideas  and  techniques.  These  explorations  have 
led  to  new  results  both  in  estimation  theory  and  in  these  other  fields,  e.g.,  to  new 
algorithms  for  triangular  and  QR  factorization  of  structured  matrices,  new  techniques 
for  root  location  and  stability  testing,  and  new  recursions  for  orthogonal  polynomials 
on  the  unit  circle  and  the  real  line  as  well  as  on  other  curves. 

For  several  years,  the  guiding  principle  in  these  studies  has  been  the  concept  of 
generalized  displacement  structure  (Lev-Ari  and  Kailath  (1986)),  which  genualizcd 
and  subsumed  our  earlier  work  on  Toeplitz-oriented  displacement  structure  (Kailath, 
Kung  and  Morf,  (1979);  see  also  Lev-Ari  and  Kailath  (1984)).  A  related  notion  of 


displacement  structure  has  also  emerged  from  recent  work  of  Heinig  and  Rost 
(1984,1987)  of  East  Germany.  While  they  are  aware  of  our  work,  and  make  some 
attempts  to  relate  to  it,  their  approach  and  methodology  are  significantly  different  from 
ours.  In  particular,  they  focus  only  on  the  problem  of  inversion  of  structured  matrices 
via  algebraic  methods,  while  our  work  has  primarily  addressed  triangular  factorization 
of  such  matrices,  and  our  approach  is  based  on  a  generating  function  characterization 
of  matrices.  The  triangular  factorization  problem  is  in  many  senses  more  fundamental 
than  inversion,  and  has  more  consequences  for  signal  processing,  linear  algebra, 
operator  theory  and  other  fields.  In  fact,  recently  we  were  able  to  reduce  the  inversion 
problem  for  structured  matrices  to  the  factorization  of  certain  block-matrices  with 
structured  blocks  (see  Chun  and  Kailath  (1989)).  This  result  also  confirms  and  clarifies 
an  earlier  observation  [see  Lev-Ari  and  Kailath  (1984)]  on  the  relation  between 
efficient  inversion  and  efficient  factorization  of  structured  matrices:  only  some  of  the 
structured  matrices  that  admit  an  efficient  factorization  procedure  can  also  be 
efficiently  inverted. 

The  generating  function  approach,  which  was  introduced  in  the  Ph.D.  research  of 
Lev-Ari  (1983),  also  suggests  a  natural  system-theoretic  interpretation  of  the  theory, 
which  allows  a  study  of  various  problems  in  m  theory,  such  as  minimal 
realization,  Pade  approximation,  control  design,  and  a  variety  of  root  distribution 


(stability)  problems  for  polynomials.  This  approach  also  unveils  the  great  flexibility  in 
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the  computational  details  of  the  factorization  procedure  for  structured  matrices.  [In 
contrast,  the  approach  taken  bv  Heinig  and  Rost  leads  to  a  single  procedure  for  (the 


inversion  of)  every  particular  type  of  structured  matrices].  In  particular,  we  were  able  ,  ifmy 
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Routh-Hurwitz  algorithms,  as  particular  instances  of  our  factorization  procedure. 
Moreover,  by  exploiting  the  flexibility  in  our  prototype  procedure  we  obtained  new 
alternatives  to  these  classical  algorithms,  with  reduced  computational  requirements  (see 
Lev-Ari,  Bistritz  and  Kailath  (1989)). 

Perhaps  the  most  prominent  system-theoretic  aspect  of  our  efficient  factorization 
techniques  is  that  they  can  be  interpreted  as  recursive  identification  procedures  for 
certain  lossless  cascade  models.  For  instance,  the  classical  Schur  algorithm  is  also  a 
procedure  for  step  by  step  identification  and  ‘peeling-  of  the  layers  of  a  transmission¬ 
line  with  a  piecewise  constant  characteristic  impedance  (Bruckstein  and  Kailath 
(1987)].  Such  layered  models  have  been  used  tor  ,;ie  some  time  in  oil  exploration 
and  in  marine  seismography.  They  involve  two  scalar  signals  propagating  in  opposite 
directions;  consequently,  the  characteristics  of  the  model  can  be  captured  by  a  single 
scalar  input-scattering  function.  Schur’s  original  formulation  of  his  algorithm  was,  in 
fact,  in  terms  of  this  scattering  function. 

We  have  recently  begun  to  extend  our  methods  to  multichannel  cascade  models, 
which  involve  involve  multiple  signals  propagating  both  in  the  forward  and  in  the 
opposite  (backward)  direction.  Since  such  models  are  represented  by  matrix  scattering 
functions,  it  would  seem  that  the  corresponding  layer-peeling  procedures  need  to  be 
rederived  in  matrix  (or  block)  form.  This  is  certainly  possible  (see,  e.g.,  Delsarte, 
Genin  and  Kamp  (1979)),  but  results  in  the  introduction  of  computation-intensive 
matrix  operations,  such  as  the  evaluation  of  matrix  square  roots.  In  contrast  to  this 
approach,  we  have  succeeded  in  obtaining  layer-peeling  procedures  that  involve  only 
elementary  (2x2)  circular  and  hyperbolic  rotations,  and  therefore  require  only  scalar 


-  4  - 


computations.  Furthermore,  our  procedures  can  be  implemented  in  pipelined  parallel 
processing  hardware  (such  as  systolic  arrays).  Consequently,  the  throughput  of  such 
implementations  is  independent  of  the  number  of  channels  (i.e.,  the  number  of  forward 
and  backward  signals  flowing  through  the  model). 

More  details  on  this  problem,  and  some  of  the  results  mentioned  above,  are  given 
in  the  following  sections. 

2.  FACTORIZATION  OF  STRUCTURED  MATRICES 

Our  early  work  on  factorization  and  inversion  of  Toeplitz  and  close-to-Toeplitz 
matrices  led  us  to  the  observation  that  for  certain  matrices  the  displacement  matrix 

VZR  :=  R  -  2RZ-  ,  Z  =  [  5,  J+I  ],"J=0 

has  low  rank.  Notice  that  Z  has  unity  elements  on  the  first  subdiagonal  and  zeros 
elsewhere.  Consequently,  the  displacement  matrix  VZR  is  the  difference  between  the 
matrix  R  and  the  matrix  ZRZ*  obtained  by  displacing  R  one  step  along  the  main 
diagonal.  In  particular,  the  displacement  rank  (i.e.,  the  rank  of  VZR)  is  2  for  both 
Toeplitz  matrices  and  inverses  of  Toeplitz  matrices.  We  have  shown  in  previous  work 
(largely  supported  by  AFOSR)  that  the  displacement  concept  is  a  key  tool  for 
developing  fast  algorithms  of  many  kinds,  including  factorization  and  inversion  of 
Toeplitz  and  near-Toeplitz  matrices,  as  well  as  fast  (generalized  Levinson  and  Schur) 
algorithms  for  solving  linear  systems  with  such  coefficient  matrices.  Not  surprisingly, 
these  results  led  naturally  to  cascade  orthogonal  structures  for  the  prediction  of 
nonstationary  processes  (Lev-Ari  and  Kailath  (1984)).  We  have  also  found  that  the 
same  concept  is  tightly  connected  to  the  more  general  problem  of  cascade  filter 
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synthesis  in  network  theory  and  digital  filtering  as  well  as  to  a  variety  of  inverse 
scattering  problems  (some  references  are  Rao  and  Kailath,  (1984,  1985),  Bruckstein 
and  Kailath  (1987)). 

Later  we  extended  the  concept  of  displacement  structure  to  a  very  broad  family  of 
structured  matrices,  including  Hankel  matrices  and  their  inverses,  sums  of  Toeplitz  and 
Hankel  matrices  and  several  others  (Lev-Ari  and  Kailath,  (1986)).  The  generalized 
displacement  of  a  matrix  R,  is  defined  as  d(Z,Z)R  where 

d(A,B)R  :=  £  *kj  A*  R  (»*)'  ,  (1) 

k,l=  o 

and  the  asterisk  (* )  denotes  Hermitian  transpose  (complex  conjugate  for  scalars). 

The  concept  of  displacement  structure  and  its  properties  are  more  conveniently 
described  in  terms  of  generating  functions.  The  generating  function  of  a  matrix  R  is 
a  power  series  in  two  complex  variables,  viz., 

R(z,w)  :=  [  1  z  z2  ...]R[1  w  w2  ...]*  (2) 


The  displacement  d(Z,Z)R  of  a  matrix  has  the  generating  function  d(z  ,w)R  (z  ,w), 
where 

d (z ,\v )  =  X  dk.i  z*(w*y  (3a) 

k,l  =0 

Thus  the  generating  function  of  a  Hermitian  matrix  with  a  displacement  structure  has 
the  form 


=  G  (z)  J  G*(w) 
d  (z  ,w ) 


(3b) 


where  J  is  any  constant  nonsingular  Hermitian  matrix.  The  triple  [d(z  ,w),G  (z)J  ) 
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is  called  a  generator  of  R{z,\v),  since  it  uniquely  determines  the  generating  function 
R(z,w). 

We  have  extended  our  previous  work  (Lev-Ari  and  Kailath  (1984))  on  efficient 
factorization  of  matrices  with  displacement  structure  to  accommodate  the  generalized 
displacement  d( Z,Z)R,  and  we  have  shown  (Lev-Ari  and  Kailath  (1986),  Lev-Ari 
(1989)),  that  efficient  factorization  of  R  is  possible  if,  and  only  if,  there  exist  power 
series  4>(z  ),\j/(z )  (with  arbitrary  radii  of  convergence)  such  that 

d  (z  ,w )  =  <J)(z  )<t>*(w )  -  y(z  )\j/*(w )  .  (4) 

To  obtain  the  factorization  of  R  one  has  to  propagate  the  recursion  (with 
Gq(z)  =  G(z)) 


(z-Q)Gi+l(z)  =  Cl(z)0,(z)  i  =  0,1,2,  .  .  .  (5a) 


where 


0;  (z ) 


I  - 


d{z  ,X  i) 


t/(z,C,¥(C,,x() 


JMi  ^  Ui 


(5b) 


M,  :=  Gi%i)RrHti&)Gi&i)  =  M;  ,  (5c) 

U t  is  any  constant  matrix  such  that 

Ut  J  U  *  =  J  ,  (5d) 

and  X(-  is  any  complex  constant  such  that 

d  (x,-  ,x,- )  =  0  .  (5e) 

The  standard  choice  of  the  extraction  points  {Q  },  i.e.,  ^  =0,  produces  triangular 
factorizations;  other  choices  can  be  useful  in  root-location  and  filter  synthesis 
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procedures  (see,  e.g.,  Deprettere  and  Dewilde  (1980),  Vaidyanathan  and  Mitra  (1984)). 

This  algorithm  requires  0(n2)  computations  to  factor  a  structured  n  xn 
matrix  R  in  the  form  R  =  LDL*,  in  contrast  to  the  conventional  LDL *  algorithm 
which  requires  0(n 3)  operations  to  factor  an  arbitrary  n  xn  matrix.  The  z'-th 
element  of  the  diagonal  matrix  D  and  the  z-th  column  of  the  lower  triangular  matrix 
L  are  determined  by  the  coefficients  of  the  power  series  expansion  of  G,(z). 


The  displacement  representation  (3)  also  leads  to  the  following  fundamental 
result:  for  every  d(z,w)  of  the  form  (4)  and  for  every  (finite)  matrix  R, 

In  [d( Z,Z)R}  =  In  {d(Z,Z)R_tl} 


(6) 


where  the  reversed  matrix  R  ^  is  obtained  by  transposing  R  1  with  respect  to  the 
anti-diagonal,  namely  the  reversed  matrix  is  defined  as 

:=  I  Qr  I  (7a) 

where  Q  can  be  any  square  matrix,  the  superscript  T  denotes  the  conventional 

(non-Hermitian)  transpose,  and  I  is  the  anti-diagonal  unity  matrix,  viz.. 


I 


(7b) 


The  fundamental  result  (6)  implies  that  R  and  R-^  have  the  same  displacement 
structure.  It  does  not  tell  us,  however,  how  to  relate  the  parametrizations  of  generators 
of  R  and  of  R-^.  Such  relations  are  known  for  Toeplitz  matrices  and  give  rise  to 
the  Levinson  algorithm  and  the  Gohberg-Semencul  formula.  They  have  been  extended 
by  Lev-Ari  and  Kailath  (1984)  to  close-to-Toeplitz  matrices,  i.e.,  to  matrices  with 
displacement  structure  involving  d(z,w)  =  1  -  zw*. 
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We  have  recently  shown  that  our  factorization  procedure  can  be  extended  to 
matrices  for  which  the  generalized  displacement  R  -  FRF*  has  low  rank,  where  F 
can  be  an  arbitrary  lower  triangular  matrix  (Chun  and  Kailath  (1989)).  [When  F  is  a 
lower-triangular  Toeplitz  matrix  the  new  procedure  reduces  to  the  one  reported  in 
(Lev-Ari  and  Kailath  (1986))].  Moreover,  we  have  shown  how  to  embed  the  problem 
of  inverting  a  structured  matrix  in  a  factorization  problem  of  a  double  size.  More 
specifically,  given  { G^/  )  such  that 

R  -  FRF*  =  GJ  G*  (8a) 

we  introduce  the  block  matrix 


M  := 


R  I 

II  O. 


(8b) 


Next  we  look  for  matrices  F,G  that  satisfy  the  generalized  displacement  equation 


r  ■> 

R  I 

O 

R  I 

F  0 

* 

G 

G 

- 

= 

J 

.1  o. 

lo  FJ 

.1  o. 

.0  F. 

G. 

.G. 

(8c) 


which  is  an  embedding  of  (8a).  The  factorization  procedure  can  now  be  used  to  obtain 
a  displacement  representation  of  the  Schur  complement  of  the  (2,2)  block  in  M,  viz., 


0  0 

F  O' 

0  0 

F  O' 

* 

r  > 

0 

O' 

- 

_ 

= 

J 

0  -R-' 

- 

lO  Fj 

0  -R-] 

^  - 

.0  F 

II. 

IL 

In  other  words,  we  obtain  a  generator  of  R_1.  A  comprehensive  account  of  this 
approach  for  the  case  F  =  Z  is  available  in  [Chun  and  Kailath  (1989)],  where  we 
have  also  shown  how  to  modify  this  method  to  obtain  extensions  of  the  Gohberg- 
Semencul  formula.  The  Ph.D.  dissertation  of  Chun  (1989)  describes  several  other 
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applications  of  suitable  embedding  procedures. 


3.  BEZOUTIANS  AND  EFFICIENT  ZERO-LOCATION 

One  application  of  the  notion  of  generalized  displacement  structure  is  the 
construction  of  Bezoutian  matrices.  Originally,  such  matrices  were  associated  with 
stability  (and  zero-location)  tests  with  respect  to  the  unit  circle  and  the  imaginary  axis. 
We  have  extended  this  notion  to  structured  matrices  whose  generating  functions  have 
the  form 

S(z.w,  =  /><«»•(■*> -K'fcHXW  (,4a) 

d(z.w) 

where 


d(z,w)  :=  a  +  pz  +(Pvv)*  +  5ziv*  ,  (14b) 

and  the  sharp  (#)  denotes  a  suitably  defined  polynomial  transformation.  The 
remarkable  property  of  such  matrices  is  that  most  of  their  elements  vanish  except  the 
elements  in  the  leading  nxn  principal  submatrix ,  where  n  :=  deg  p(z).  This 
submatrix,  which  we  denote  by  B,  has  full  rank  (i.e.,  rank  B  =  n)  if,  and  only  if, 
p{z)  and  pn(z)  are  coprime.  Furthermore,  the  inertia  of  the  Bezoutian  matrix  B 
(i.e.,  the  number  of  its  positive,  zero  and  negative  eigenvalues)  serves  to  locate  the 
zeros  of  the  polynomial  p{z)  with  respect  to  the  following  partition  of  the  complex 
plane, 

Q+  -  [z  ;  d(z,z)  >  0) 

:=  {z  ■  d  (z  ,z)  —  0)  (15) 

Q_  :=  {z  ;  d(z  ,z )  <  0) 


More  specifically  the  inertia  of  B  indicates  how  many  zeros  are  shared  by  p(z)  and 
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pn(z)  and  how  many  of  the  remaining  zeros  are  in  Q+  and  in  Q_  .  Our  fast 
factorization  procedure  makes  it  possible  to  determine  the  inertia  of  a  Bezoutian  matrix 
in  0(n2)  opera^  .ms,  starting  with  p(z)  and  ptt(z),  and  without  explicit  evaluation 
of  the  elements  of  B.  For  Bezoutians  on  the  imaginary  axis  and  the  unit  circle  our 
formulation  leads  (among  other  possibilities)  to  the  Routh-Hurwitz  and  Schur-Cohn 
tests,  and  serves  to  delimit  the  family  of  O  (n2)  polynomial  zero-location  procedures 
(Lev-Ari,  Bistritz  and  Kailath  (1989)). 

4.  RECURSIVE  LAYER  PEELING 

The  fundamental  factorization  procedure  (5)  for  structured  matrices,  viz., 

(z  -Ci)Gi+1(z)  =  Gi(r)0i(z) 

is,  at  the  same  time,  also  a  layer-peeling  procedure.  Starting  with  G0(z),  which  w'e 
can  interpret  as  boundary  data  for  a  layered  medium,  we  identify  an  elementary  layer 
with  chain-scattering  matrix  0o(z),  then  “peel”  it  off  to  obtain  Gflz),  the 
boundary  data  for  the  rest  of  the  medium  (with  the  first  layer  removed),  and  repeat  the 
same  procedure  again  and  again.  Such  layer-peeling  recursions  have  been  used  in 
cascade  filter  synthesis  (see,  e.g.,  Dewilde,  Vieira  and  Kailath,  (1978);  Vaidyanathan 
and  Mitra  (1984)),  in  inverse  scattering  (Bruckstein  and  Kailath  (1987)),  zero-location 
(Lev-Ari,  Bistritz  and  Kailath  (1987)),  and  model-order  reduction  (see,  e.g.,  Genin  and 
Kung  (1981)). 

The  classical  work  of  Schur  (1917)  forms  the  basis  for  much  of  the  subsequent 
work  on  layer  peeling  procedures.  Schur’s  algorithm  associates  a  sequence  of  so- 
called  reflection  coefficients,  all  with  magnitude  bounded  by  unity,  with  every’  passive 
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scattering  function,  i.e.,  a  function  /  (z)  that  is  analytic  and  bounded  by  1  in  the 
unit  disc.  In  particular,  if  /(z)  is  an  all-pass  function ,  which  means  that 
| y  (e-7 e)  |  =1  for  all  9,  then  Schur’s  algorithm  produces  a  finite  sequence  of 
reflection  coefficients  {&,-  ;  0  <i  <n),  where  \kn  |  =1  and  |  kt  |  <1  for 
0  <  i  <  n  -  I. 

Another  property  of  the  algorithm  is  that  starting  with  a  passive  scattering 
function  it  generates  a  sequence  of  such  functions.  This  is  the  essence  of  the  layer 
peeling  method:  a  single  step  of  the  Schur  algorithm  applied  to  a  passive  medium 
leaves  a  medium  with  the  same  property,  which  makes  it  possible  to  apply  the  same 
step  again  and  again.  A  single  step  of  the  Schur  algorithm  corresponds  to  the  removal 
(or  peeling)  of  an  elementary  lossless  two-port.  Thus,  the  algorithm  produces  a 
discrete  transmission-line  model,  whose  input  scattering  function  is  /(z)  (Fig.  1): 


Figure  1.  Transmission-line  model  associated  with  the  Schur  algorithm. 

In  addition  to  the  recursive  characterization  of  passivity  via  the  constraint  on  the 
magnitude  of  the  reflection  coefficients,  Schur  also  introduced  an  operator-norm 
characterization  of  passivity:  he  proved  that  for  any  function  /(z)  that  is  analytic  in 
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the  unit  disc  we  can  construct  an  infinite  lower-triangular  Toeplitz  matrix  whose  first 
column  consists  of  the  cocffiripnts  of  the  power  series  expansion  of  /  (z ),  viz., 

/o 

/.  / 0  0 

W)  =  fi  f  i  •  06a) 

such  that 

sup  |/(z)|  =  H  L(/ )  II 2  <  1  (16b) 

lz  I  <  1 

where  \\A  ||2  denotes  the  conventional  (spectral)  norm  of  a  matrix  A,  i.e., 

|| Ax  ||2 

\\A  ||  2  ■-  sup  — (16c) 

*  lUII 2 

and  |1  jc  ||  2  denotes  the  Euclidean  (/2)  norm  of  a  vector  x. 

Schur’s  algorithm  was  later  extended  by  Cohn  (1922)  to  functions  with  poles  in 
the  unit  disc,  but  only  to  rational  “all-pass”  functions,  i.e.,  to  functions  / (?)  of  the 
form 

/(Z)  =  X£^  (17a) 

P(z) 

where  p#  (z)  denotes  the  conjugate  reversal  operation,  viz., 

pn(z)  =  zdes  P^[p(l/z*)Y  (17b) 

The  now  well-known  Schur-Cohn  test  associates  with  each  such  function+  a  finite 
sequence  of  reflection  coefficients,  some  of  which  have  magnitudes  larger  than  1. 
Moreover,  it  has  been  shown  (e.g.,  using  the  properties  of  Bezoutians  on  the  unit  disc) 
that  the  number  of  poles  of  /  (z)  =  pn  (z)lp{z)  inside  the  unit-disc  equals  the  number 

+  Assuming  p(z)  has  no  «;ros  at  Z  =  0,  and  applying  the  algorithm  to  / (z  )/ X. 
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of  singular  values  of  the  matrix  L (f )  that  are  larger  than  |  X  |  or,  equivalently,  the 
number  of  negative  eigenvalues  of  the  following  finite  rank  matrix, 

R:=  |  X 1 2I  -  L(  /  )L*(  /  )  .  (18) 

Rational  allpass  functions  of  a  given  degree  k  are  members  in  the  family  H£° 

of  all  functions  with  k  poles  (or  less)  inside  the  unit  circle,  and  whose  magnitude  is 

bounded  on  the  unit  circle,  i.e.,  sup  |/(z)|  <  °°.  It  turns  out  that  the  Schur-Cohn 

lz  I  =  I 

algorithm  does  not  map  the  family  H£°  into  itself,  except  when  k=0.  This  means 
that  this  algorithm  does  not  admit  the  same  layer-peeling  interpretation  as  the  standard 
Schur  algorithm.  Nevertheless,  we  have  found  that  it  is  possible  to  modify  the  Schur 
algorithm  in  such  a  way  that  the  resulting  recursion  indeed  maps  the  family  H£°  into 
itself  and,  therefore,  admits  the  same  layer-peeling  interpretation  as  the  classical  Schur 
algorithm.  Moreover,  the  layers  involve  only  elementary  (2x2)  orthogonal  and 
hyperbolic  rotations  (Ackner,  Lev-Ari  and  Kailath  (1989)). 

Furthermore,  our  modified  recursion  applies  to  every  function  f(z)eH£°,  and 
not  just  to  allpass  functions.  Each  layer  in  the  resulting  transmission-line  model  has  an 
sign  or  ‘polarity’  (Fig.  2).  While  a  p-type  layer  (i.e.,  one  of  positive  polarity)  maps 
//t“  into  itself,  a  n-type  layer  maps  H£°  into  namely,  it  reduces  by  one  the 

number  of  poles  within  the  unit  disc.  Therefore,  the  number  of  n-lype  layers  in  the 
transmission-line  model  that  is  generated  by  our  modified  algorithm  equals  the  number 
of  poles  that  the  function  /(z)  has  within  the  unit  disc.  This  is  the  key  idea  in  the 
construction  of  efficient  procedures  for  zero-location  and  model-order  reduction,  and  in 
the  solution  of  several  interpolation  and  extension  problems  arising  in  various  control 
applications  (Ackner,  Lev-Ari  and  Kailath  (1989)). 


L 
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p-type  layer  (kn  =  /„( 0))  n-type  layer  ( kn  =  l//*(0)) 

Figure  2.  Layers  for  the  modified  Schur  algorithm 

Recursive  layer  peeling  for  a  medium  with  multiple  inputs  and  outputs  involves  a 
generalization  of  the  Schur  algorithm  to  matrix-valued  analytic  functions.  One  version 
of  the  matrix  Schur  algorithm  (Delsarte,  Genin  and  Kamp  (1979))  requires  hyperbolic 
matrix  rotations,  which  are  computationally  expensive  since  they  involve  the  finding  of 
the  square  root  of  (positive  definite)  matrices.  Moreover,  this  approach  cannot  be 
generalized  to  matrix  functions  with  elements  in  //t",  because  now  it  would  involve 
square  roots  of  indefinite  matrices. 

An  alternative  approach  to  the  matrix  analytic  case  was  taken  by  Dewilde  and 
Dym  (1981).  They  extract  simpler  layers  than  in  the  method  of  Delsarte,  Genin  and 
Kamp  and,  as  a  result,  there  is  no  need  to  compute  square  roots  of  matrices.  This 
computational  procedure,  which  is  known  as  a  tangential  Schur  algorithm,  was  the 
starting  point  for  our  research  on  layer-peeling  methods  for  MIMO  (multiple- 
input/multiple-output)  systems. 


We  have  shown  that  the  procedure  of  Dewilde  and  Dym  can  be  transformed  into 
an  equivalent  form  that  involves  only  elementary  (2x2)  rotations.  In  fact,  for  pxq 
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matrix  scattering  functions,  the  peeling  of  each  layer  in  our  version  of  the  algorithm  is 
implemented  by  a  sequence  of  q- 1  elementary  orthogonal  rotations,  a  single 
hyperbolic  rotation,  and  a  single  block-delay  (Fig.  3).  These  operations  can  be  easily 
implemented  in  either  software  or  hardware. 


Figure  3.  Single  layer  of  the  multichannel  Schur  algorithm  (p  =  3  ,  q  =4) 

In  addition  to  requiring  significantly  fewer  computations  than  the  tangential  Schur 
algorithm,  our  version  also  serves  to  clarify  the  relationship  between  the  MIMO  case 
and  the  better-known  scalar  or  SISO  (single-input/single-output)  case.  The  MIMO 
procedure  differs  from  the  scalar  one  only  by  the  presence  of  elementary  orthogonal 
rotations.  Thus  both  procedures  share  the  same  core,  which  consists  of  elementary 
hyperbolic  rotations  and  delay  elements.  Consequently,  we  were  able  to  show  that  our 
version  of  the  tangential  Schur  algorithm  can  be  modified  to  accommodate  scattering 
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functions  with  poles  within  the  unit  circle,  and  that  this  modification  affects  only  the 
‘core’  part  of  the  procedure  (i.e.,  it  is  independent  of  the  number  of  inputs  or  outputs). 
This  means  that  each  layer  still  has  a  ‘sign’  or  polarity,  as  in  the  SISO  case,  and  that 
the  number  of  n-type  layers  coincides  with  the  number  of  Smith-McMillan  poles  of 
the  given  matrix  scattering  function.  This  work  is  part  of  the  dissertation  research  of  R. 
Ackner,  to  be  completed  by  June  1990  or  earlier. 
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